subroutine jagp_local_energy_integral_der_aa_aaaa_f( &
       ni, & ! (in)    number of occupied orbitals
       na, & ! (in)    number of unoccupied orbitals
       RT, & ! (in)    unoccupied row matrix times inverse of occ pairing matrix
      RTt, & ! (in)    transpose of RT
      Jsi, & ! (in)    single index occ jastrow intermediate
      Jsa, & ! (in)    single index unocc jastrow intermediate
     Jdia, & ! (in)    double index occ-unocc jastrow intermediate
     Jdai, & ! (in)    double index unocc-occ jastrow intermediate
     Jdij, & ! (in)    double index occ-occ jastrow intermediate
     Jdab, & ! (in)    double index unocc-unocc jastrow intermediate
      Tia, & ! (in)    occ-unocc slice of the one-electron integrals
    Viajb, & ! (in)    occ-unocc-occ-unocc slice of the two-electron integrals
        E, & ! (out)   energy contribution
      Wia, & ! (out)   adjoint for occ-unocc slice of the one-electron integrals
    Wiajb)   ! (out)   adjoint for occ-unocc-occ-unocc slice of the two-electron integrals

implicit none

integer,      intent(in)    ::    ni
integer,      intent(in)    ::    na
real(kind=8), intent(in)    ::    RT(na,ni)
real(kind=8), intent(in)    ::    RTt(ni,na)
real(kind=8), intent(in)    ::    Jsi(ni)
real(kind=8), intent(in)    ::    Jsa(na)
real(kind=8), intent(in)    ::    Jdia(ni,na)
real(kind=8), intent(in)    ::    Jdai(na,ni)
real(kind=8), intent(in)    ::    Jdij(ni,ni)
real(kind=8), intent(in)    ::    Jdab(na,na)
real(kind=8), intent(in)    ::    Tia(ni,na)
real(kind=8), intent(in)    ::    Viajb(ni,na,ni,na)
real(kind=8), intent(out)   ::    E
real(kind=8), intent(out)   ::    Wia(ni,na)
real(kind=8), intent(out)   ::    Wiajb(ni,na,ni,na)

integer :: i, a, j, b

real(kind=8) :: tempQ

real(kind=8) :: yjb, yj, yaj, yab, ya, yia, yij, yib, yi, yiajb

! initialize energy to zero
E = 0.00d+00

! one electron integrals
! begin a-loop
do a = 1,na

  ! begin i-loop
  do i = 1,ni

    ! compute intermediates
    yia = Jsa(a) * Jdia(i,a)
    yi  = yia * Jsi(i)

    ! get jastrow-modified integral
    yiajb = Tia(i,a) * yi

    ! get the contribution to the local energy apart from the jastrow-modified integral
    tempQ = RTt(i,a)

    ! add local energy contribution
    E = E + tempQ * yiajb

    ! get integral adjoint derivative
    Wia(i,a) = yi * tempQ

  ! end i-loop
  enddo

! end a-loop
enddo


! two electron integrals
! begin b-loop
do b = 1,na

  ! begin j-loop
  do j = 1,ni

    ! compute j-loop intermediates
    yjb = Jsa(b) * Jdia(j,b)
    yj  = yjb * Jsi(j)

    ! begin a-loop
    do a = 1,na

      ! compute a-loop intermediates
      yab = yj  * Jdab(a,b)
      yaj = yab * Jdai(a,j)
      ya  = yaj * Jsa(a)

      ! begin i-loop
      do i = 1,ni

        ! compute i-loop intermediates
        yib = ya  * Jdia(i,b)
        yij = yib * Jdij(i,j)
        yia = yij * Jdia(i,a)
        yi  = yia * Jsi(i)

        ! get jastrow-modified tei element
        yiajb = Viajb(i,a,j,b) * yi

        ! get the portion of the energy contribution other than the jastrow-modified integral
        tempQ = 0.50d+00 * ( RTt(i,a) * RTt(j,b) - RT(a,j) * RTt(i,b) )

        ! add energy contribution
        E = E + yiajb * tempQ

        ! get integral adjoint derivative
        Wiajb(i,a,j,b) = yi * tempQ

      ! end i-loop
      enddo

    ! end a-loop
    enddo

  ! end j-loop
  enddo

! end b-loop
enddo

end subroutine jagp_local_energy_integral_der_aa_aaaa_f




subroutine jagp_local_energy_integral_der_aabb_f( &
       ni, & ! (in)    number of occupied orbitals
       na, & ! (in)    number of unoccupied orbitals
      RTt, & ! (in)    transpose of contraction of unoccupied rows with inverse matrix
       TC, & ! (in)    contraction of unoccupied columns with inverse matrix
      RTC, & ! (in)    contraction of RT with unoccupied columns
     IPMt, & ! (in)    transpose of the inverse pairing matrix
      UPM, & ! (in)    the unoccupied pairing matrix
      Jsi, & ! (in)    single index jastrow intermediate for i
      Jsa, & ! (in)    single index jastrow intermediate for a
      Jsj, & ! (in)    single index jastrow intermediate for j
      Jsb, & ! (in)    single index jastrow intermediate for b
     Jdia, & ! (in)    double index jastrow intermediate for i and a
     Jdij, & ! (in)    double index jastrow intermediate for i and j
     Jdib, & ! (in)    double index jastrow intermediate for i and b
     Jdaj, & ! (in)    double index jastrow intermediate for a and j
     Jdab, & ! (in)    double index jastrow intermediate for a and b
     Jdjb, & ! (in)    double index jastrow intermediate for j and b
    Viajb, & ! (in)    two electron integrals in 1122 order
        E, & ! (out)   energy contribution
    Wiajb)   ! (out)   adjoint for the current two electron integral slice

implicit none

integer,      intent(in)    ::   ni
integer,      intent(in)    ::   na
real(kind=8), intent(in)    ::   RTt(ni,na)
real(kind=8), intent(in)    ::   TC(ni,na)
real(kind=8), intent(in)    ::   RTC(na,na)
real(kind=8), intent(in)    ::   IPMt(ni,ni)
real(kind=8), intent(in)    ::   UPM(na,na)
real(kind=8), intent(in)    ::   Jsi(ni)
real(kind=8), intent(in)    ::   Jsa(na)
real(kind=8), intent(in)    ::   Jsj(ni)
real(kind=8), intent(in)    ::   Jsb(na)
real(kind=8), intent(in)    ::   Jdia(ni,na)
real(kind=8), intent(in)    ::   Jdij(ni,ni)
real(kind=8), intent(in)    ::   Jdib(ni,na)
real(kind=8), intent(in)    ::   Jdaj(na,ni)
real(kind=8), intent(in)    ::   Jdab(na,na)
real(kind=8), intent(in)    ::   Jdjb(ni,na)
real(kind=8), intent(in)    ::   Viajb(ni,na,ni,na)
real(kind=8), intent(out)   ::   E
real(kind=8), intent(out)   ::   Wiajb(ni,na,ni,na)

integer :: i, a, j, b

real(kind=8) :: tempQ, tempARTC

real(kind=8) :: yjb, yj, yaj, yab, ya, yia, yij, yib, yi, yiajb

! initialize energy to zero
E = 0.00d+00

! begin b-loop
do b = 1,na

  ! begin j-loop
  do j = 1,ni

    ! compute j-loop intermediates
    yjb = Jsb(b) * Jdjb(j,b)
    yj  = yjb * Jsj(j)

    ! begin a-loop
    do a = 1,na

      ! compute a-loop intermediates
      yab = yj  * Jdab(a,b)
      yaj = yab * Jdaj(a,j)
      ya  = yaj * Jsa(a)
      tempARTC = UPM(a,b) - RTC(a,b)

      ! begin i-loop
      do i = 1,ni

        ! compute i-loop intermediates
        yib = ya  * Jdib(i,b)
        yij = yib * Jdij(i,j)
        yia = yij * Jdia(i,a)
        yi  = yia * Jsi(i)

        ! get jastrow-modified integral
        yiajb = Viajb(i,a,j,b) * yi

        ! remember the rest of the local energy term
        tempQ = RTt(i,a) * TC(j,b) + IPMt(i,j) * tempARTC

        ! add energy contribution
        E = E + yiajb * tempQ

        ! get integral adjoint derivative
        Wiajb(i,a,j,b) = yi * tempQ

      ! end i-loop
      enddo

    ! end a-loop
    enddo

  ! end j-loop
  enddo

! end b-loop
enddo

end subroutine jagp_local_energy_integral_der_aabb_f
